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Abstract The Q-modeling strategy is followed to derive a new subgrid parame- 
terization of the turbulent stress tensor in large-eddy simulation (LES). 
The LES-Q modeling yields an explicitly filtered subgrid parameteri- 
zation which contains the filtered nonlinear gradient model as well as 
a model which represents 'Leray-regularization'. The LES-q model is 
compared with similarity and eddy- viscosity models that also use the dy- 
namic procedure. Numerical simulations of a turbulent mixing layer are 
performed using both a second order, and a fourth order accurate finite 
volume discretization. The Leray model emerges as the most accurate, 
robust and computationally efficient among the three LES-Q subgrid 
parameterizations for the turbulent mixing layer. The evolution of the 
resolved kinetic energy is analyzed and the various subgrid-model con- 
tributions to it are identified. By comparing LES-Q at different subgrid 
resolutions, an impression of finite volume discretization error dynamics 
is obtained. 
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1. Introduction 

Accurate modeling and simulation of turbulent flow is a topic of in- 
tense ongoing research. The approaches to this problem area can be 
distinguished, e.g., by the amount of detail that is intended to be in- 
cluded in the physical and numerical description. Simulation strategies 
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that aim to calculate the full, unsteady solution of the governing Navier- 
Stokes equations are known as direct numerical simulations (DNS). The 
DNS approach does not involve any modeling or approximation except 
its numerical nature and in principle it can provide solutions that possess 
all dynamically relevant flow features [1, 2]. In turbulent flow, these fea- 
tures range from large, geometry dependent scales to very much smaller 
dissipative length-scales. While accurate in principle, the DNS approach 
is severely restricted by limitations in spatial and temporal resolution, 
even with modern computational capabilities, because of the tendency 
of fluid flow to cascade its energy to smaller and smaller scales. 

This situation summons alternative, restricted simulation approaches 
to the turbulent flow problem that are aimed at capturing the primary 
features of the flow above a certain length-scale only. A prominent exam- 
ple of this is the large-eddy simulation (LES) strategy [3]. Rather than 
aiming for a precise and complete numerical treatment of all features 
that play a role in the evolution of the flow, an element of turbulence 
modeling is involved in LES [4]. In the filtering approach to LES, this 
modeling element is introduced by applying a spatially localized filter op- 
eration to the Navier-Stokes equations [5] . This introduces a smoothing 
of the flow features and a corresponding reduction in the flow complexity 
[6]. One commonly adopts spatial convolution filters which effectively 
remove the small-scale flow features that fall below an externally intro- 
duced length-scale A, referred to as the filter- width. This smoothing can 
significantly reduce the requirements on the resolution and, thus, allow 
LES to be performed for much more realistic situations than DNS, e.g., 
at higher Reynolds number, within the same computational capabilities 
[7]. This constitutes the main virtue of LES. 

The LES approach is conceptually different from the Reynolds Aver- 
aged Navier-Stokes (or, RANS) approach, which is based on statistical 
arguments and exact ensemble averages that raise the classic turbulence 
closure problem. When the spatially localized smoothing operation in 
LES is applied to the nonlinear convective terms in the Navier-Stokes 
equations, this also gives rise to a closure problem that needs to be re- 
solved. Thus, the LES approach must face its own turbulence closure 
problem: How to model the effects of the filtered-out scales in terms of 
the remaining resolved fields? 

In the absence of a comprehensive theory of turbulence, empirical 
knowledge about subgrid-scale modeling is essential but still incomplete. 
Since in LES only the dynamical effects of the smaller scales need to be 
represented, the modeling is supposed to be simpler and more straight- 
forward, compared to the setting encountered in statistical modeling 
such as in RANS. To guide the construction of suitable models we ad- 
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vocate the use of constraints based on rigorous properties of the LES 
modeling problem such as realizability conditions [8] and algebraic iden- 
tities [5, 6]. A thoughtful overview of these constraints is given in [9]. 

In this paper, we follow the a-modeling approach to the LES closure 
problem. The a-modeling approach is based on the Lagrangian-averaged 
Navier-Stokes— a equations (LANS— a, or NS— a) described below. The 
LANS— a approach eliminates some of the heuristic elements that would 
otherwise be involved in the modeling. The original LANS— a theory also 
involves an elliptic operator inversion in defining its stress tensor. When 
we apply filtering in defining the LANS— a stress tensor, instead of the 
operator inversion in the original theory, we call it LES-a . 

Background and references for LANS— a, or NS— a equations. 

The inviscid LANS— a equations (called Euler— a, in the absence of 
viscosity) were introduced through a variational formulation in [10], 
[11] as a generalization to 3D of the integrable inviscid ID Camassa- 
Holm equation discovered in [12]. A connection between turbulence and 
the solutions of the viscous 3D Camassa-Holm, or Navier-Stokes-alpha 
(NS— a) equations was identified, when viscosity was introduced in [13]- 
[15]. Specifically, the steady analytical solution of the NS— a equations 
was found to compare successfully with experimental and numerical data 
for mean velocity and Reynolds stresses for turbulent flows in pipes and 
channels over a wide range of Reynolds numbers. These comparisons 
suggested the NS— a equations could be used as a closure model for the 
mean effects of subgrid excitations. Numerical tests further substanti- 
ating this intuition were performed and reported in [16]. 

An alternative more "physical" derivation for the inviscid NS— a equa- 
tions (Euler— a), was introduced in [17] (see also [14]). This alterna- 
tive derivation was based on substituting in Hamilton's principle the 
decomposition of the Lagrangian fluid-parcel trajectory into its mean 
and fluctuating components at linear order in the fluctuation amplitude. 
This was followed by making the Taylor hypothesis for frozen-in turbu- 
lence and averaging at constant Lagrangian coordinate, before taking 
variations. Hence, the descriptive name Lagrangian-averaged Navier- 
Stokes— a equations (LANS— a) was given for the viscous version of this 
model. A variant of this approach was also elaborated in [18] but this 
resulted in a second-grade fluid model, instead of the viscous LANS— a 
equations, because the choice of dissipation made in [18] differed from the 
Navier-Stokes dissipation chosen in [13]-[17]. The geometry and anal- 
ysis of the inviscid Euler— a equations was presented in [19], [20]. The 
analysis of global existence and well-posedness for the viscous LANS— a 
was given for periodic domains in [21] and was modified for bounded 
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domains in [22]. For more information and a guide to the previous lit- 
erature specifically about the NS— a model, see paper [23]. The latter 
paper also discusses connections to standard concepts and scaling laws in 
turbulence modeling, including relationships of the NS— a model to large 
eddy simulation (LES) models that are pursued farther in the present 
paper. Related results interpreting the NS— a model as an extension of 
scale similarity LES models of turbulence are also reported in [24]. A 
numerical comparison of LANS— a model results with LES models for 
the late stages of decaying homogeneous turbulence is discussed in [25]. 
Vortex interactions in the early stages of 3D turbulence decay are stud- 
ied numerically with LANS— a and compared with both DNS and the 
Smagorinsky eddy viscosity approach in [26]. 

Three contributions in the present approach. Stated most sim- 
ply, the LANS— a approach may be interpreted as a closure model for the 
turbulent stress tensor that is derived from Kelvin's circulation theorem, 
using a smoothed transport velocity, as discussed in [23], [24], [27]. A 
new development within this approach is introduced here that gives rise 
to an explicitly filtered similarity- type model [28] for the turbulent stress 
tensor, composed of three different contributions. The first contribution 
is a filtered version of the nonlinear gradient model. The unfiltered 
version of this model is also known as the 'Clark' model [29, 30], the 
'gradient' model [31] or the 'tensor-diffusivity' model [32]. The second 
contribution, when combined with the filtered nonlinear gradient model, 
represents the so-called 'Leray regularization' of Navier-Stokes dynamics 
[33] . Finally, a new third contribution emerges from the derivation which 
completes the full LES-a model and endows it with its own Kelvin's cir- 
culation theorem. 

To investigate the physical and numerical properties of the resulting 
three-part LANS— a subgrid parameterization, we consider a turbulent 
mixing layer [34] . This flow is well documented in literature and provides 
a realistic canonical flow problem [35] suitable for testing and compar- 
ison with predictions arising from more traditional subgrid model de- 
velopments [7]. In particular, we consider similarity, and eddy- viscosity 
modeling, combined with the dynamic procedure based on Germano's 
identity [5, 6, 36], to compare with LES-a . In addition to the full 
LES-a model, in our comparisons we also consider the two models that 
are contained in it, i.e., the filtered nonlinear gradient model and the 
Leray model. We will refer to all three as LES-a models. For all these 
models, the explicit filtering stage is essential. Without this filtering op- 
eration in the definition of the models, a finite time instability is observed 
to arise in the simulations. The basis for this instability can be traced 
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back to the presence of antidiffusion in the nonlinear gradient contri- 
bution. We sketch an analysis of the one-dimensional Burgers equation, 
following [31], to illustrate this instability and show, through simulation, 
that increasing the subgrid resolution further enhances this instability 
Analyzing the resolved kinetic energy dynamics reveals that this insta- 
bility is associated with an excessive contribution to back-scatter. 

The 'nonlinearly dispersive' filtered models that arise in LES-a are 
reminiscent of similarity LES models [24]. The LES-a model separates 
the resolved kinetic energy (RKE) of the flow into the sum of two con- 
tributions: namely, the energies due to motions at scales that are either 
greater, or less than an externally determined length-scale (a). The two 
contributions are modeled by 

RKE = RKE& + RKEW (1) 

As we shall describe later in reviewing the LES-a strategy the kinetic 
energy RKE^ of turbulent motions at scales less than a is modeled 
by a term proportional to the rate of dissipation of the kinetic energy 
RKE^ at scales greater than a. (The time-scale in the proportionality 
constant is the viscous diffusion time a 2 /2v.) 

A key aspect of the LES-a dynamics is the exchange, or conversion, of 
kinetic energy between RKE^ and RKE^\ We focus on the contri- 
butions to the dynamics of the resolved kinetic energy RKE» at scales 
greater than a that arise from the different terms in the LES-a models. 
The filtered nonlinear gradient model, the Leray model and the full 
LES-a model all contribute to the reduction of the RKE^ in the lam- 
inar stages of the flow. This corresponds to forward scatter of RKE^ 
into RKE^ outweighing backward scatter. In the developing turbulent 
flow regime, the resolved kinetic energy RKE& of the full LES-a model 
may decrease too slowly compared to DNS, and for some settings of 
the (numerical) parameters can even become reactive in nature, thereby 
back-scattering too much kinetic energy from RKE^ into RKE^ . In 
contrast, the contribution of the Leray model to the RKE^ dynamics 
remains forward in nature and appears to settle around some negative, 
nonzero value in the turbulent regime. All the LES-a models show con- 
tributions to both forward and backward scatter of RKE. It is observed 
that in the full LES-a model two of the three terms almost cancel in the 
evolution of resolved kinetic energy RKE^\ This cancellation nearly 
reduces the full LES-a model to the filtered nonlinear gradient model. 

The mixing layer simulations indicate that the Leray subgrid model 
provides more accurate predictions compared to both the filtered non- 
linear gradient and the full LES-a model. This is based on comparisons 
that include mean flow quantities, fluctuating flow properties and the 
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energy spectrum. In addition, the Leray LES-a model appears more 
robust with respect to changes in numerical parameters. Predictions 
based on this model compare quite favorably with those obtained using 
dynamic (mixed) models and filtered DNS results. The Leray model 
combines this feature with a strongly reduced computational cost and is 
favored for this reason, as well. In addition, a number of classic math- 
ematical properties (e.g., existence and uniquenes of strong solutions) 
can be proven rigorously for fluid flows that are modeled with Leray's 
regularization. These can be used to guide further developments of this 
model such as extensions to more complex flows at higher Reynolds num- 
ber. This is a topic of current research and will be published elsewhere 

[37]. ' , . , , ' , 

Apart from the problem of modeling the subgrid-scale stresses, any 
actual realization of LES is inherently endowed with (strongly) inter- 
acting errors arising from the required use of marginal numerical res- 
olution [38, 39, 40, 41]. The accuracy of the predictions depends on 
the numerical method and subgrid resolution one uses. We consider in 
some detail numerical contamination of a 'nonlinear gradient fluid' and 
a 'Leray fluid,' which are defined as the hypothetical fluids governed by 
the corresponding subgrid model. In this analysis we are consequently 
not concerned with how accurately the modeled equations represent fil- 
tered DNS results. Rather, we focus on the numerical contamination of 
the predictions. For this purpose we compare two finite volume spatial 
discretization methods, one at second order, and the other at fourth 
order accuracy. 

The subgrid modeling and the spatial discretization of the equations 
give rise to a computational dynamical system whose properties are in- 
tended to simulate those of the filtered Navier-Stokes equations. The 
success of this simulation depends of course on the properties of the 
model, as well as of the spatial discretization method and the subgrid 
resolution. The model properties are particularly important in view of 
the marginal subgrid resolution used in present-day LES. We consider 
the role of the numerical method at various resolutions and various ratios 
of the filter-width A compared to the grid-spacing h. Let A be a fixed 
constant. In cases of large ratios A//i > 1 one approximates the grid- 
independent LES solution corresponding to the given value of A, and 
the accuracy of its predictions will be limited by the quality of the as- 
sumed subgrid model. At the other extreme, one may assume A/h to be 
rather small and numerical effects can constitute a large source of error. 
Through a systematic variation of the ratio A/h at constant A we can 
identify the contributions of the numerical method at coarse resolutions. 
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This will give an impression of how the computational dynamical system 
is affected by variations in the resolution and the numerical method. 

The organization of this chapter is as follows. In section 2 we intro- 
duce the large-eddy simulation problem and identify the closure problem 
and some of its properties. The treatment of this closure problem using 
the a-framework is sketched, together with more conventional subgrid 
parameterization that involves the introduction of similarity, and eddy- 
viscosity modeling. Finally, we analyze the instabilities associated with 
the use of the unfiltered nonlinear gradient model. In section 3 we in- 
troduce the numerical methods used and consider the simulation of a 
turbulent mixing layer. Some direct and large-eddy simulation results 
will be shown. In section 4 we focus on the LES-a models and consider 
the dynamics of the resolved kinetic energy in each of the three cases. 
This comparison provides a framework for understanding how the differ- 
ent LES-a subgrid models function. We proceed with an assessment of 
the numerical error dynamics at relatively coarse subgrid resolutions. A 
summary and concluding remarks for the chapter are given in section 5. 

2. Large-eddy simulation and a-modeling 

This section sketches the traditional approach to large-eddy simu- 
lation, which arises from direct spatial filtering of the Navier-Stokes 
equations (section 2.1). The algebraic and analytic properties of the 
LES modeling problem will be discussed first. The LES closure problem 
will then be considered in the a-framework of turbulent flow, derived 
via Kelvin's circulation theorem for a smoothed, spatially filtered trans- 
port velocity (section 2.2). The closure of the filtered fluid flow problem 
achieved this way will be compared with the more traditional methods of 
similarity, and eddy- viscosity modeling for LES. The latter is introduced 
in section 2.3 together with the dynamic procedure based on Germano's 
identity. We also sketch a stability analysis of the one-dimensional fil- 
tered Burgers equation involving the nonlinear gradient subgrid model 
that illustrates the instabilities associated with this model (section 2.4). 

2.1. Spatially filtered fluid dynamics 

We consider the incompressible flow problem in d spatial dimensions. 
The Cartesian velocity fields Uj (i = 1, . . . , d) and the normalized pres- 
sure field p constitute the complete solution. The velocity field is con- 
sidered to be solenoidal and the evolution of the solution is described by 
the Navier-Stokes equations. These are conservation laws for mass and 
momentum, respectively, that can be written in the absence of forcing 
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as 

djUj = (2) 
d t Ui + dj{uiUj) + dip - j^djjUi = (3) 

where dt and dj denote, respectively, partial differential operators in time 
t and Cartesian coordinate Xj, j = 1, . . . , d. The quantity Re = u r l r /v r 
is the Reynolds number based on reference velocity (u r ), reference length 
(l r ) and reference kinematic viscosity (v r ), which were selected to non- 
dimensionalize the governing equations. Repeated indices are summed 
over their range, except where otherwise noted. 

Equations (2 - 3) model incompressible flow in all its spatial and 
temporal details. In deriving approximate equations that are specialized 
to capture the generic large-scale flow features only, one applies a spatial 
filter operation L : u — > u to (2 - 3). For simplicity, we restrict to linear 
convolution filters: 

u(x, t) = L(u) (x, t) = f°° G(x - £) u(£, t) d£=(G* u) (x, t) (4) 

in which the filter-kernel G is normalized, i.e., L(c) = c for any constant 
solution u = c. We assume that the filter-kernel G is localized as a 
function of x— £ and a filter- width A can be assigned to it. Typical filters 
which are commonly considered in LES are the top-hat, the Gaussian 
and the spectral cut-off filter. Here, we restrict ourselves to the top-hat 
filter which has a filter-kernel given by 

G(z) = { A ~ 3 if ^ <Ai/2 (5) 
^ ' \ otherwise ^ ' 

where Aj denotes the filter-width in the X{ direction and the total filter- 
width A is specified by 

A 3 = A!A 2 A 3 (6) 

Apart from the filter-kernel in physical space, the Fourier-transform of 
G(z), denoted by H(k), is important, e.g., for the interpretation of the 
effect of the filter-operation on signals which are composed of various 
length-scales. The Fourier-transform of the top-hat filter is given by (no 
sum in Ajfej) 

If we consider a general Fourier-representation of a solution u(x, t), 

n(x,t) = ^ Ck (t)e ik - x (8) 

k 
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the filtered solution can directly be written as 

Tl(x,t) = ^(i7(k) Ck (t)) e tkx (9) 

k 

We notice that each Fourier-coefficient c^(t) is attenuated by a fac- 
tor -ff(k). The normalization condition of the filter-operation implies 
H(Q) = 1. For small values of |Ajfej| one infers from a Taylor expansion 

H(h) = 1 - (1/24) ((^A^ 2 + (k 2 A 2 f + (k 3 A 3 ) 2 ) +... (10) 

which shows the small attenuation of flow features which are consider- 
ably larger than the filter-width A, i.e., |Ajfej| <C 1 for % = 1, 2, 3. As |k| 
increases H(k) becomes smaller while oscillating as a function of Aj/uj. 
Consequently, the coefficients Ck(t) are strongly reduced as |Ajfcj| 3> 1 
and the small scale features in the solution are effectively taken out by 
the filter operation. Similarly, the Gaussian filter can be shown to have 
the same expansion for small |Ajfej| and reduces to zero monotonously 
as |Ajfcj| becomes large. 

The filter operation L is a convolution integral. Hence, it is a linear 
operation that commutates with partial derivatives [42, 43]. This prop- 
erty facilitates the application of the filter to the governing equations (2 
- 3). A straightforward application of such filters leads to 

djUj = (11) 

d t Ui + djiuiUj) + dip - -^-djjUi = -djTij (12) 

tie 

where we introduced the turbulent stress tensor 

Tij = UjUj — UiUj (13) 

We observe that the filtered solution {Hi, p} represents an incompressible 
flow (djUj = 0). The same differential operator as in (3) acts on {ui, p} 
and due to the filtering a non-zero right-hand side has arisen which 
contains the divergence of the turbulent stress tensor r^. This latter 
term is the so-called subgrid term, and expressing it in terms of the 
filtered velocity and its derivatives constitutes the closure problem in 
large-eddy simulation. 

The LES modeling problem as expressed above has a number of im- 
portant, rigorous properties which may serve as guidelines for specifying 
appropriate subgrid-models for t^. In particular, we will briefly review 
realizability conditions, algebraic identities and transformation proper- 
ties. Adhering to these basic features of limits some of the heuristic 
elements in the subgrid-modeling. 
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Realizability. It is well known that the Reynolds stress u^u'j in 
RANS is positive semi-definite [44, 45] and the following inequalities 
hold [46] 

m > for ie {1,2,3} (no sum) (14) 

\nj\ < y/TuTjj for i,je {1,2,3} (no sum) (15) 
det(ry) > (16) 



If the filtering approach is followed, in general r^- ^ u^u'j and, therefore, 
it is relevant to know the conditions under which is positive semi- 
definite. Following Vreman et al. [8], it can be proved that in LES 
is positive semi-definite if and only if the filter kernel G(x, £) is positive 
for all x and £. If we assume G > 0, the expression 

(f,g)= ! G(x,Z)f{Z)g(S)dt (17) 
Jn 

defines an inner product and we can rewrite the turbulent stress tensor 
as: 

^■(x) = f G(x,0(«i(0 - = K,vf) (18) 

with vf(£) = Uj(^) — Uj(x). In this way the tensor forms a 3 x 3 
Grammian matrix of inner products. Such a matrix is always positive 
semi-definite and consequently r%j satisfies the realizability conditions. 
The reverse statement can likewise be established, showing that the con- 
dition G > is both necessary and sufficient. 

One prefers the turbulent stress tensor in LES to be realizable for 
a number of reasons. For example, if is realizable, the generalized 
turbulent kinetic energy k — Ta/1 is a positive quantity. This quantity 
is required to be positive in subgrid models which involve the fc-equation 
[47]. Several further benefits of realizability and positive filters can be 
identified [8] ; here we restrict to adding that the kinetic energy of u is 
bounded by that of u for positive filter- kernels: 

If If 

: 'u\ 2 dx < - / |ti| 2 (ix (19) 



2 Jn ~ 2 Jn 

Requiring realizability places some restrictions on subgrid models. For 
example, if G > models for r should be realizable. Consider, e.g., an 
eddy-viscosity model given by 

2 

mj = -VeOij + -^k5ij (20) 
In order for this model to be realizable, a lower bound for k in terms of 



the eddy- viscosity v e arises, i.e., k > ^VSaVe where a = ^ij^ij and Oij 
is the rate of strain tensor given by = diUj + djUi. 
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Algebraic identities. The introduction of the product operator 
S(u{,Uj) = UiUj allows to write the turbulent stress tensor as [6]: 

T-j = uiu]-UiUj = L(S(ui,Uj)) - S(L(ui),L(uj)) = [L,S](ui,Uj) (21) 

in terms of the central commutator [L, S] of the filter L and the product 
operator S. This commutator shares a number of properties with the 
Poisson-bracket in classical mechanics. Leibniz' rule of Poisson-brackets 
is in the context of LES known as Germano's identity [5] 

[£ 1 £ 2 ,S} = [£ 1 ,S}£ 2 + £ 1 [£ 2 ,S] (22) 

This can also be written as 

t £i£j = t Ci C 2 + C 1 t C2 (23) 

and expresses the relation between the turbulent stress tensor corre- 
sponding to different filter-levels. In these identities, C\ and C 2 denote 
any two filter operators and t k = [K, S] . The first term on the right- 
hand side of (23) is interpreted as the 'resolved' term which in an LES 
can be evaluated without further approximation. The other two terms 
require modeling of r at the corresponding filter-levels. 
Similarly, Jacobi's identity holds for S, £\ and L 2 : 

[C u [£ 2 , S}\ + [£ 2 , [S, £i]] = - [S, [A, C 2 }} (24) 

The expressions in (22) and (24) provide relations between the turbulent 
stress tensor corresponding to different filters and these can be used to 
dynamically model t l . The success of models incorporating Germano's 
identity (22) is by now well established in applications for many different 
flows. In the traditional formulation one selects C\ = H and C 2 = L 
where H is the so called test-filter. In this case one can specify Germano's 
identity [5] as 

t hl {u) = t h {L{u)) + H (r L (u)) (25) 

The first term on the right hand side involves the operator t h acting on 
the resolved LES field L(u) and during an LES this is known explicitly. 
The remaining terms need to be replaced by a model. In the dynamic 
modeling [36] the next step is to assume a base-model m K correspond- 
ing to filter-level K and optimize any coefficients in it, e.g., in a least 
squares sense [48]. The operator formulation can be extended to include 
approximate inversion defined by £~ 1 (L(x k )) = x k for < k < N [49]. 
Dynamic inverse models have been applied in mixing layers [50]. 
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Transformation properties. The turbulent stress tensor r^- can be 
shown to be invariant with respect to Galilean transformations. This 
property also holds for the divergence, i.e., djTij, referred to as the sub- 
grid scale force. Hence, the filtered Navier-Stokes dynamics is Galilean 
invariant. Suitable subgrid models should at least maintain the Galilean 
invariance of the divergence of the model, i.e., djiriij. In fact, most 
subgrid models are represented by tensors which are Galilean invariant. 
Examples of non-symmetric tensor models have been reported in [24] for 
which, however, djiriij was verified to be Galilean invariant. 

Likewise, it is of interest to consider a transformation of the subgrid 
scale stress tensor to a frame of reference rotating with a uniform angular 
velocity. The full subgrid scale stress tensor transforms in such a way 
that the subgrid scale force is the same in an inertial and in a rotating 
frame. Horiuti [51] has recently analyzed several subgrid scale models 
and showed that some of them do not satisfy this condition. This is an 
example of how transformation properties of the exact turbulent stress 
tensor can be used to guide propositions for subgrid modeling. 

After closing the filtered equations (11-12) by a subgrid model stress 
tensor rriij we arrive at the modeled filtered dynamics described in the 
absence of forcing by 

d jVj = (26) 

dm + dj(viVj) + diP - —djjVi = -djiriij (27) 

tie 

whose solution is denoted as {vi, P}. Ideally, if m^- and the numerical 
treatment were correct and had no undesirable effects on the dynamics, 
one might expect i>j = U{. In view of possible sensitive dependence of 
an actual solution, e.g., on the initial condition, one should not expect 
instantaneous and point-wise equality of Vi and Ui but rather one should 
expect statistical properties of the filtered and modeled solution to be 
equal. Assessing the extent to which the properties of {v j, P} and {ui,p} 
are correlated allows an evaluation of the quality of the subgrid model, 
the dynamic effects arising from the numerical method and the interac- 
tions between modeling and numerics. In what follows, we will use the 
notation {vi,P} to distinguish the solution of the subgrid model from 
the filtered solution {lli,p}. 

2.2. Subgrid model derived from Kelvin's 
theorem 

The LES-a modeling scheme we shall use here is based on the well- 
known viscous Camassa-Holm equations, or LANS— a model. This mod- 
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eling strategy imposes a "cost" in resolved kinetic energy (RKE) for 
creation of smaller and smaller excitations below a certain, externally 
specified length scale, denoted by a. This cost in converting RKE^ to 
RKE^ implies a nonlinear modification of the Navier-Stokes equations 
which is reactive, or dispersive, in nature instead of being diffusive, as is 
more common in present-day LES modeling. The modification appears 
in the nonlinear convection term and can be rewritten in terms of a sub- 
grid model for the turbulent stress tensor. In the LANS— a model, the 
processes of nonlinear conversion of RKE» to RKE«> and sweeping 
of the smaller scales by the larger ones are still included in the modeled 
dynamics. We will sketch the LES-a approach in this subsection and 
extract the subgrid models used in this study. For more details and 
applications of this approach, see [13]-[17], [23]-[26]. 

It is well known that the Navier-Stokes equations satisfy Kelvin's cir- 
culation theorem, i.e., 

<£ Uj dxj = I d kk Uj dxj (28) 

dt / 7 (u) 7 7 (u) Re 

Here 7(11) represents a fluid loop that moves with the Eulerian fluid 
velocity u(x, t). The basic equations in the LES-a modeling may be 
introduced by modifying the velocity field by which the fluid loop is 
transported. The governing LES-a equations will provide the smoothed 
solution {vj, P} and we specify the equations for v through the Kelvin- 
filtered circulation theorem. Namely, we integrate an approximately 'de- 
filtered' velocity w around a loop 7(v) that moves with the regularized 
spatially filtered fluid velocity v, cf. [23], [24], [27] 

Jti iv) WjdXj = i iv) T e dkkW > dx > w 

Hence, the basic transport properties of the LES-a model arise from 
filtering the 'loop-velocity' to obtain v, then approximately defiltering v 
to obtain the velocity w in the Kelvin integrand. As we shall show, this 
approach will yield the model stress tensor needed to complete the 
filtering approach outlined in section 2.1. Direct calculation of the time 
derivative in this modified circulation theorem yields the Kelvin-filtered 
Navier-Stokes equations, 

d t Wi + VjdjWi + w k div k + diP - -^-djjWi = , djvj = (30) 

He 

where we introduce the scalar function P in removing the loop integral. 
The relation between the 'defiltered' velocity components Wi and the 
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LES-a velocity components Vi of the Kelvin loop needs to be specified 
separately. The Helmholtz defiltering operation was introduced in [10], 
[11] for this purpose: 

Wi = Vi- a 2 djjVi = (1 - a 2 djj)Vi = H a (vi) (31) 

where 7i a denotes the Helmholtz operator. We recall that all explicit 
filter operations L with a non-zero second moment, have a Taylor ex- 
pansion whose leading order terms are of the same form as (31). Conse- 
quently, we infer that the leading order relation between a and A follows 
as a 2 = A 2 /24 for the top-hat and the Gaussian filter. We will use this 
as the definition of a in the sequel. 

The LES-a equations can be rearranged into a form similar to the ba- 
sic LES equations (27), by splitting off a subgrid model for the turbulent 
stress tensor. For the Helmholtz defiltering, we obtain from (30): 

d t Vi + dj(viVj) + diP + djmfj - —djjVi = , djVj = (32) 

after absorbing gradient terms into the redefined pressure P. Thus, we 
arrive at the following parameterization for the turbulent stress tensor 

Haimfj) = a 2 (d k Vi d k vj + d k Vi djV k - diV k djV^j (33) 

In the evaluation of the LES-a dynamics in the above formulation, an 
inversion of the Helmholtz operator Ti a is required. The 'exponential' 
(or 'Yukawa') filter [52] is the exact, explicit filter which inverts TL a . 
Thus, an inversion of 7i a corresponds to applying the exponential filter 
to the right-hand side of (33) in order to find mfj. However, since the 
Taylor expansion of the exponential filter is identical at quadratic order 
to that of the top-hat and the Gaussian filters, we will approximate the 
inverse of 7i a by an application of the explicit top-hat filter, for reasons 
of computational efficiency. Moreover, in actual simulations the numeri- 
cal realization of the exponential filter is only approximate and can just 
as well be replaced by the numerical top-hat filter. This issue of (ap- 
proximately) inverting the Helmholtz operator will be studied separately 
and published elsewhere [37] . 

The full LES-a subgrid model mfj has three distinct contributions. 
The first term on the right-hand side is readily recognized as the nonlin- 
ear gradient model which we will denote by A; L j . This term is closely re- 
lated to the similarity model proposed by Bardina [28], as will be shown 
in the next subsection. The second term will be denoted by Bij and 
combined with the first term, corresponds to the Leray regularization of 
the convective terms in the Navier-Stokes equations. This regularization 
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arises if the familiar contribution UjdjUi in the Navier-Stokes equations 
is replaced by UjdjVi in the smoothed description. The third term will 
be denoted by Cij. Further details of the derivation and mathematical 
properties of the LES-a model will be published elsewhere [37]. We can 
explicitly write the stress tensor for the LES-a model as 



The explicit filter, represented by the overbar in this expression, is real- 
ized by the numerical top-hat filter in this study. It does not necessarily 
have to coincide with the LES-filter. While the LES-filter specifies the 
relation between the Navier-Stokes solution Ui and the LES-a solution 
Vi, the explicit LES-a filter is used to approximate TL~ l . We will con- 
sider the effects associated with variations in the filter-width A of the 
explicit LES-a filter with filter- width A/ A = k. Typical values that will 
be considered are k = 1 and k = 2. 

In the next subsection we will describe some familiar subgrid mod- 
els used in LES which are based on the similarity and eddy-viscosity 
concepts. 

2.3. Similarity modeling and eddy-viscosity 



We distinguish two main contributions in present-day traditional sub- 
grid modeling of the turbulent stress tensor, i.e., dissipative and simi- 
larity subgrid models. In this subsection we briefly describe these two 
basic approaches, as well as subgrid models that consist of combinations 
of an eddy- viscosity and a similarity part, so-called mixed models. The 
relative importance of the two components in such mixed models is ob- 
tained by using the dynamic procedure which is based upon Germano's 
identity (25). This mixed approach effectively regularizes and stabilizes 
similarity models. 

As a result of the filtering, flow features of length-scales (much) smaller 
than the filter-width A are considerably attenuated. This implies that 
the natural molecular dissipation arising from the viscous fluxes, is 
strongly reduced, compared to the unfiltered flow-problem. In order 
to compensate for this, dissipative subgrid-models have been introduced 
to model the turbulent stress tensor. The prime example of such eddy- 
viscosity models is the Smagorinsky model [2, 53]: 




(34) 



regularization 



m, 



! 3 = -(C s A)V(v)Mv) with Kv)| 2 = I^-(v)^(v) (35) 



16 



where is the strain rate, introduced above (crjj = diVj + djVi). This 
model adds only little computational overhead. The major short-coming 
of the Smagorinsky model is its excessive dissipation in laminar regions 
with mean shear, because is large in such regions [36]. Furthermore, 
the correlation between the Smagorinsky model and the actual turbulent 
stress is quite low (reported to be ~ 0.3 in several flows). 

In trying to compensate for these short-comings of the Smagorinsky 
model, a second main branch of subgrid models emerges from the simi- 
larity concept [28]. Using the commutator notation, the turbulent stress 
tensor can be expressed as Ty(u) = [L, S](ui,Uj). In terms of this short- 
hand notation, the basic similarity model can be written as 

mfj = [L,S](vi,Vj) (36) 

i.e., directly following the definition of the turbulent stress tensor, but 
expressed in terms of the available modeled LES velocity field. Gen- 
eralizations of this similarity model arise by replacing Vi in (36) by an 
approximately defiltered field % = C{vi) where C(L(u)) « u, i.e., C 
approximates the 'inverse' of the filter L [49]. In detail, a generalized 
similarity model arises from mP = [L, S] using the approximate 

inversion. This approach is also known as the deconvolution model [54] 
and is reminiscent to the subgrid estimation model [55]. The correlation 
with Tjj is much better with correlation coefficients reported in the range 
0.6 to 0.9 in several flows. The low level of dissipation associated with 
these models renders them quite sensitive to the spatial resolution. At 
relatively coarse resolutions, the low level of dissipation can give rise to 
instability of the simulations. Moreover, these models add significantly 
to the required computational effort. At suitable resolution, however, 
the predictions arising from generalized similarity models are quite ac- 
curate. 

An interesting subgrid model which follows the similarity approach to 
some degree and avoids the costly additional filter-operations is the non- 
linear gradient model, mentioned earlier. This model can be derived from 
the Bardina scale-similarity model by using Taylor expansions of the fil- 
tered velocity. One may arrive at r^- = ^ Y^k ^{d k Ui){d k Uj) + C(A 4 ). 
The first term on the right-hand side is referred to as the 'nonlinear 
gradient model' or tensor-diffusivity model: 

rnl D = ^ A Y J tf k {d k v i ){d k v j ) (37) 

k 

Since this model is part of all three different LES-a models identified in 
the previous subsection, we will analyze the dynamics and the instabili- 
ties arising from this model in some more detail in the next subsection. 
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The three subgrid models, i.e., (35), (36) and (37) constitute well- 
known examples in LES-literature, which represent basic dissipative and 
reactive, or dispersive, properties of subgrid models for the turbulent 
stress tensor. These basic similarity and eddy-viscosity models can be 
combined in mixed models using the dynamic procedure, which provides 
a way of combining the two basic components of a mixed model without 
introducing additional external ad hoc parameters. 

We consider simple mixed models based on eddy-viscosity and similar- 
ity. In these models the eddy-viscosity component reflects local turbu- 
lence activities and the local value of the eddy-viscosity adapts itself to 
the instantaneous flow. The dynamic procedure starts from Germano's 
identity (25). A common way to write Germano's identity is: 

Tij = Rij (38) 

where 

= Ujuj — HiUj (39) 

Rij = (uiUj)-UiUj (40) 

Here, in addition to the basic LES-filter (•) of width A a so-called 'test'- 
filter (•) of width A is introduced. Usually, this test-filter is wider than 
the LES-filter and the combined filter (•) has a width that follows from 

A = A 2 +A . This relation is exact for the composition of two Gaussian 
filters and can be shown to be 'optimal' for other filters such as the top- 
hat filter [7]. The only external parameter that needs to be specified in 
the dynamic procedure is the ratio A/ A which is commonly set equal 
to two. The terms at the left-hand side of the Germano identity (38) 
are the turbulent stress tensor on the 'combined' filter level (Ty) and 
the turbulent stress tensor, filtered with the test-filter (r^), respectively. 
Finally, Rij represents the resolved stress tensor which can be explicitly 
calculated using the modeled LES fields. 

The general procedure for obtaining 'locally' optimal model param- 
eters in a mixed formulation starts by assuming a basic model rriij to 
approximate the turbulent stress tensor r^, and a corresponding model 
Mij for T^. We consider rriij to be of 'mixed' type, i.e., 

rriij = aij + cbij (41) 

where and bij are basic models. These basic models involve opera- 
tions on v only; = ajj(v), bij = hj(v). Furthermore, in standard 
mixed models, c is a scalar coefficient-field which is to be determined. 
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The model Mjj is represented as: 

Mij = Aij + CBij (42) 

where = ajj(v), B^ = 6y(v). It is essential in this formulation 
that the coefficient C corresponding to the composed filter-level is well 
approximated by the coefficient c; i.e., we assume C ~ c. Insertion in 
Germano's identity yields + m^- = or in more detail, 



Aij + cijj ) + c 



-^27 "I - bij 



Rij (43) 



where we have used the approximation cbij ~ c6y. Introducing the 
short-hand notation Aij = A; L j + &ij, Bij = B, L j + bij, the coefficient c 
is required to obey cBij = Rij — Aij. This relation should hold for all 
tensor-components, which of course is not possible for a scalar coefficient 
field c. To resolve this situation we introduce an averaging operator (/) 
and define the 'Germano-residual' by 

e{c) = ( l -{{R i j-A i j)-cBij} 2 ) (44) 

From this we obtain an optimality condition for c from e'(c) = and we 
can solve the local coefficient as 

_ ((Rjj -Aij)Bij) 

~ (Bi&j) {4b) 

where we assumed (cfg) w c(fg). The averaging operator (/) is usually 
defined in terms of an integration over homogeneous directions of the 
flow-domain. In the case of the mixing layer, considered here, the aver- 
aging over the homogeneous streamwise and spanwise direction results 
in a dynamic coefficient c which is a function of the normal coordinate 
X2 and time t. In more complex flow-domains, averaging over homoge- 
neous directions may no longer be possible. Taking a running-average 
over time t is then a viable alternative, as was recently established, e.g., 
for flow in a spatially developing mixing layer [35]. 

As an example we consider the Smagorinsky model as the base model. 
The corresponding models on the two filter-levels can be written as 



m: 



-C d A>(v)|<ry(v) ; M% = -C d A \a(v)\aij(v) (46) 



The 'optimal' Cd follows from = (Rij Bij) / ' (BijBij) . In order to pre- 
vent numerical instability caused by negative values of Cd, the model 
coefficient Cd is artificially set to zero at locations where the proce- 
dure would return negative values. Sometimes, in developing flows, it 
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is beneficial to also introduce a 'ceiling'-value for This value should 
be chosen such that once the flow is well-developed in time the actual 
limitation arising from the ceiling- value is no longer restrictive [35]. 

The dynamic mixed model employs the sum of Bardina's similarity 
and Smagorinsky eddy- viscosity model as the base model, i.e., 

= [L, S](v t , Vj ) - C d A 2 |5(v)|^(v) (47) 
Likewise, a mixed nonlinear gradient model can be introduced by 

m% G = 1a 2 (cM(^-) " C d A 2 |5(v)|^(v) (48) 

The dynamic procedure has been used in a number of different flows. 
Compared to predictions using only the constitutive base models, the dy- 
namic procedure generally enhances the accuracy and robustness. More- 
over, it responds to the developing flow in such a way that the eddy- 
viscosity is strongly reduced in laminar regions and near solid walls [56] . 
This avoids specific modeling of transitional regions and near-wall phe- 
nomena, provided the resolution is sufficient. At even coarser resolution 
one may have to resort to specific models for transition and walls. We 
will not enter into this problem. Rather, we will focus on the properties 
of the nonlinear gradient model in the next subsection. 

2.4. Analysis of instabilities of the nonlinear 
gradient model in one dimension 

From the discussion of the previous two subsections, it would appear 
that the nonlinear gradient subgrid model would be very well suited 
to parameterize the dynamic effects of the small scales in a turbulent 
flow. This model is part of the full LES-a model and it also emerges 
as a Taylor expansion of the Bardina similarity model. In this sub- 
section we will analyse the nonlinear gradient model in the context of 
the one-dimensional Burgers equation and show that this model gives 
rise to very strong instabilities. Apparently, some features appear to be 
missing in the pure nonlinear gradient model. In subsequent sections 
we will show in what way the explicit filtering and the other terms in 
the LES-a model, or dynamic eddy-viscosity regularization, alter this 
peculiar behavior of the nonlinear gradient model. 

We will analyse the nature of the instability of the pure gradient model 
for the one-dimensional Burgers equation [31]. The linear stability of a 
sinusoidal profile will be investigated. If a flow is linearly unstable then 
it is nonlinearly unstable to arbitrarily small initial disturbances. The 
linear analysis thus provides some information on the nonlinear equation. 
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The Burgers equation with gradient subgrid-model is written as: 

d t u + \d x {u 2 ) - vd 2 x u = -^ V d x (d x u) 2 + f{x) (49) 

The parameter r\ = A 2 /12. The following analysis shows that smooth 
solutions of equation (49) can be extremely sensitive to small pertur- 
bations, leading to severe instabilities. In particular, we consider the 
linear stability of a 27r-periodic, stationary solution, U (x, t) = sin(x) on 
the domain [0, 2ir] with periodic boundary conditions. The forcing func- 
tion / is determined by the requirement that U is a solution of equation 
(49). We substitute a superposition of U and a perturbation w, 

u(x,t) = U(x) +w(x,t) (50) 

into equation (49) and linearize around U, omitting higher order terms 
in w: 

dtw + (1 — rj) sin(x) d x w + (w + f]d 2 w) cos(x) = ud 2 w (51) 

We use a Fourier expansion for w written as w = «fc(i)e ifcr . After 
substitution of this series into equation (51) we obtain an infinite system 
of ordinary differential equations for the Fourier coefficients a^: 

«fe = 7; k (v k -V- l)«fc-i - k 2 ua k + \k{rjk + rj + l)a k +i (52) 



To understand the nature of the nonlinear gradient model for the 
Burgers equation, we first analyse system (52) assuming v = 0. In- 
stead of the infinite system, we consider a sequence of finite dimensional 
systems, 

(53) 



M n z n 



where z n is a vector containing the 2n + 1 Fourier coefficients a- n ...a r , 
and M n is a (2n + 1) x (2n + 1) tri-diagonal matrix: 



a. 



M n 



l n 

r n ■ 



■ h 
r 2 h 

h r 2 
h ■ 



■ r n 
In 



(54) 
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with 

h = ^Hvk - V ~ 1) ; r fc = -(fc-l)(7/fc + l) (55) 

The eigenvalues of M n determine the stability of the problem. The 
system is unstable if the maximum of the real parts of the eigenvalues 
is positive. We denote the eigenvalues of M n by Xj and introduce X max 
such that 

\X m ax\ = max|Aj| (56) 

3 

This eigenvalue problem can be shown to have the following asymptotic 
properties (for a detailed proof see [31]): 

1. if A is an eigenvalue then — A is an eigenvalue (57) 

2. \X max \ ~ r]n 2 (58) 

3. |Im(Amox)|<n-l (59) 

The first point implies that X max can be chosen such that Re(X max ) > 0. 
Hence, the combination of these three properties yields the asymptotic 
behavior of the maximum of the real parts of the eigenvalues: 

Re(A maa; ) ~ r/ra 2 (60) 

This shows that the inviscid system is linearly unstable and that the 
largest real part of the eigenvalues is asymptotically proportional to n 2 , 
where n is the number of Fourier modes taken into account. 

It should be observed that the instability is severe, since the system 
is not only unstable, but the growth rate of the instability is infinitely 
large as n — > oo. The instability is fully due to the incorporation of 
the gradient model, since all eigenvalues of the matrix M n are purely 
imaginary in case the inviscid Burgers equation without subgrid-model 
is considered (77 = 0). In numerical simulations the instability will grow 
with a finite speed, since then the number of Fourier modes is limited by 
the finite grid. Moreover, expression (60) illustrates that grid-refinement 
(with 77 kept constant), which corresponds to a larger n, will not stabilize 
the system, but rather enhance the instability. The growth rate of the 
instability of the one-dimensional problem can be expressed in terms of 
A and the grid-spacing h: rjii 2 ~ (A//1) 2 . Consequently, the instability 
is not enhanced if the ratio between A and h is kept constant. 

Finally, we will consider the more complicated case v 7^ 0. The linear 
system in equation (52) now gives rise to matrices M n which have a 
negative principal diagonal. It is known that for every fixed value of 
n there exists an eigenvalue arbitrarily close to the eigenvalue of the 
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inviscid system (X max ) if v is sufficiently small [57]. Hence for small 
values of v the viscous system for finite n is still linearly unstable. The 
matrix M n is strictly diagonally dominant if u > 77 + 1, while all rows 
except n and n+2 are already diagonally dominant if v > 77. If the matrix 
is diagonally dominant, the real parts of all eigenvalues are negative and, 
consequently, the system is stable. This indicates that stability can be 
achieved by a sufficiently large viscosity, which does not depend on n, 
but only on 77. Thus, if the gradient model is supplemented with an 
adequate eddy-viscosity the instability will be removed as is the case 
with a dynamic mixed model involving the gradient model. 

3. Numerical simulations of a turbulent mixing 
layer 

In this section we first present the numerical methods used to solve 
the DNS and LES equations (subsection 3.1). We illustrate the accuracy 
of these methods for turbulent flow in a mixing layer in subsection 3.2. 

3.1. Time-integration and spatial discretization 

The Navier-Stokes or modeled LES equations are discretized using 
the so-called method of lines. We consider the compressible formulation 
and perform simulations at a low convective Mach number which was 
shown to provide essentially incompressible flow-dynamics. The method 
of lines allows to treat the spatial and temporal discretization separately 
and gives rise to a large number of ordinary differential equations for the 
unknowns on a computational grid. 

We write the Navier-Stokes or LES equations concisely as dfU = 
where 14 denotes the state-vector containing, e.g., velocity and pressure, 
and T is the total flux, composed of the convective, the viscous, and 
possibly the subgrid fluxes. The operator T contains first and second 
order partial derivatives with respect to the spatial coordinates Xj. The 
equations are discretized on a uniform rectangular grid and the grid 
size in the Xj-direction is denoted by hj. If we adopt a specific spatial 
discretization around a grid point x^fc, the operator T{U) is approxi- 
mated in a consistent manner by an algebraic expression Fij k ({U a /3y}) 
where {U a ^} denotes the state vectors in all the grid-points, labeled by 
a, (3, 7. Usually, only neighboring grid points around (i, j, k) appear 
explicitly in Fij k , e.g., in case finite difference or finite volume discretiza- 
tions are considered. After applying the method of lines, the governing 
equations yield 

d t U ijk (t) = F ijk {{U aPl }) ; U ijk (0) = U§l (61) 
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where U-jl represents the initial condition. Hence, in order to specify 
the numerical treatment, apart from the initial and boundary conditions, 
the spatial discretization which gives rise to and the temporal inte- 
gration need to be specified. We next introduce these separately. 

The time stepping method which we adopt is an explicit four-stage 
compact-storage Runge-Kutta method. When we consider the scalar 
differential equation du/dt = f(u), this Runge-Kutta method performs 
within one time step of size St 

U U) = u (0) + PjStf^O- 1 )) (j = 1, 2, 3, 4) (62) 

with u(°) = u(t) and u(t + St) = u (4 \ With the coefficients ft = 1/4, 
ft = 1/3, ft = 1/2 and ft = 1 this yields a second-order accurate 
time integration method [58]. The time step is determined by the sta- 
bility restriction of the numerical scheme. It depends on the grid-size 
h and the eigenvalues of the flux Jacobi matrix of the numerical flux 
/. In a short-hand notation one may write St = CFL h/\X max \ where 
|A maa; | denotes the eigenvalue of the flux Jacobi matrix with maximal 
size, and CFL denotes the Courant-Friedrichs-Levy-number which de- 
pends on the specific choice of explicit time integration method. For 
the present four-stage Runge-Kutta method a maximum CFL number 
of 2.4 can be established using a Von Neumann stability analysis. In the 
actual simulations we use CFL = 1.5, which is suitable for both DNS 
and LES, irrespective of the specific subgrid model used. 

In order to specify the spatial discretization we distinguish between 
the treatment of the convective and the viscous fluxes. We will only 
specify the numerical approximation of the ft-operator; the ft and ft- 
operators are treated analogously. Subgrid-terms are discretized with 
the same method as the viscous terms. Throughout we will use a second 
order method for the viscous fluxes and both a second order, and a fourth 
order accurate method for the convective fluxes. All these methods 
are constructed from (a combination of) first order numerical derivative 
operators Dj. 

The second-order method that we consider is a finite volume method 
[59]. The discretization of the convective terms is the cell vertex trape- 
zoidal rule, which is a weighted second-order central difference. In vertex 
(i, j, k) the corresponding operator is denoted by D\ and for the approx- 
imation of d\f it is defined as 

(Dif)ij,k = (si+i,j,k ~ Si-i,j,k)/{2hi) (63) 
with s iJik = {9ij-i,k + ^9i,j,k + &j+i,fc)/4 
and g ijjik = (fij,k-i + %fi,j,k + /;j,fc+i)/4 
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The viscous terms contain second-order derivatives which are treated by 
a consecutive application of two first order numerical derivatives. This 
requires for example that the gradient of the velocity is calculated in 
centers of grid-cells. In center {i + + \,k + |) the corresponding 
discretization D 2 f has the form 

(£>2A + i J+ i, fe+ i = (64) 
with s iJ+ i )fe+ i = (fi, jt k + fi,j+i,k + fi,j,k+i + /ij+i,fe+i)/4 

The second derivative is subsequently calculated with operator D\ ; thus 
we approximate, e.g., dn(f) ijk « D 1 (D 2 (f)) ij k- 

The combination of D\ and D 2 is robust with respect to odd-even 
decoupling but it is only second order accurate. In a similar manner 
we may construct a fourth-order accurate method. The corresponding 
expression for D$f has the following form: 

{D3f)i,j,k = {~Si+2,j,k + 8sj+ij,fc ~ 8s i-i,j,fe + s i-2,j,k)/ (12/ii) (65) 
with Sij >k = (-9i,j-2,k + %,j-i,fc + Wgi,j,k + %,j+i,fc - 5 , jJ+2,fe)/16 
and = (—fij t k-2 + ^fi,j,k-i + Wfi,j,k + ^fi,j,k+i - /j,j,fc+2)/16 

This scheme is conservative, since it is a weighted central difference. 
The coefficients in the definition for gijk are chosen such that gi jk is 
a fourth order accurate approximation to fij^k and 7r-waves in the X3- 
direction give no contributions to gij k- The definition for Sij k nas t ne 
same properties with respect to the X2-direction. For convenience, we 
will refer to a combination of D3 for the convective, and D\, D2 for the 
viscous fluxes as fourth-order methods, but we remark that the formal 
spatial accuracy of the scheme is only second-order due to the treatment 
of the viscous terms. 

3.2. The turbulent mixing layer 

The flow in a temporally developing turbulent mixing layer is well 
documented in literature (e.g. [7]), and will be considered here to test 
the LES-a modeling approach. In this section we review the scenario of 
the development of the flow that is considered and sketch the type of 
predictions that can be obtained by traditional LES using the dynamic 
model. This serves as a point of reference for the next section. 

We simulate the compressible three-dimensional temporal mixing layer 
and use a convective Mach number M = 0.2 and a Reynolds number 
based on upper stream velocity and half the initial vorticity thickness 
of 50. The governing equations are solved in a cubic geometry of side 
/ = 59. Periodic boundary conditions are imposed in the streamwise 



LES-a modeling of turbulent mixing 



25 



(x±) and spanwise (x^) direction, while in the normal (x 2 ) direction 
the boundaries are free-slip walls. The initial condition is formed by 
mean profiles corresponding to constant pressure p = l/^M 2 ) where 
7 = 1.4 is the adiabatic gas constant, u\ = tanh(x2) for the streamwise 
velocity component, u 2 = 113 = and a temperature profile given by 
the Busemann-Crocco law. Superimposed on the mean profile are two- 
and three-dimensional perturbation modes obtained from linear stability 
theory. Further details may be found in [34]. 




Figure 1. Results from a DNS using 192 3 points. Contours of spanwise vorticity 
for the plane X3 = 3L/4 at t = 20, t — 40 and t = 80 from left to right. Solid and 
dotted contours indicate negative and positive vorticity respectively. The contour 
increment is 0.1. 

The DNS is conducted on a uniform grid with 192 3 cells using the 
fourth order spatial discretization method. Visualization of the DNS 
data demonstrates the roll-up of the fundamental instability and suc- 
cessive pairings (figure 1). Four rollers with mainly negative spanwise 
vorticity are observed at t = 20. After the first pairing (t = 40) the flow 
has become highly three-dimensional. Another pairing (t = 80), yields 
a single roller in which the flow exhibits a complex structure. 

The accuracy of the simulation with 192 3 cells is satisfactory as is 
inferred from coarser grid computations on 64 3 and 128 3 cells. The 
evolution of the momentum thickness 

S(t) = \ f L/2 (l-(n 1 ))((u 1 ) + l)dx 2 (66) 

4 J-L/2 

and an instantaneous velocity component at the center of the shear layer 
are shown in figure 2. The 64 3 -simulation is inadequate for the prediction 
of the local instantaneous solution, but the momentum thickness appears 
quite reasonable. 

To illustrate the effect that filtering has on a well-developed DNS 
solution, vorticity contours for A = L/16 are shown in figure 3. Com- 
paring this with the corresponding DNS results in figure 1 allows one to 
appreciate the strong smoothing effect that filtering has on the solution. 
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Figure 2. Evolution of the momentum thickness (left) and U3 at (\L, 0, \L) (right) 
obtained from simulations which do not involve any subgrid model and employ a 
sequence of grids: 64 3 (dotted), 128 a (dashed) and 192 3 (solid). 
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Figure 3. Left: Contour- lines of the 0-component of the vorticity. The effect of 
spatially filtering the DNS solution at t — 80 in figure 1 with a top-hat filter and 
filter-width A = L/16. Right: prediction of the kinetic energy with the dynamic 
eddy-viscosity model (dashed) compared with the filtered DNS results (markers) and 
a simulation on the coarse LES grid (32 3 ) without a model (solid). 

On the right in figure 3 we included the decay of the resolved turbulent 
kinetic energy, defined as 

E = \( («? + «! + «§) dx (67) 
2 Jn 

We observe that the dynamic eddy-viscosity model generates quite a 
correction of the 'no- model coarse grid simulation'. Other models were 
also considered in [7], such as Smagorinsky's model, Bardina's scale- 
similarity model and dynamic mixed models. Roughly speaking, the use 
of Bardina's model leads to flow predictions which contain somewhat 
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too many small scale features whereas the Smagorinsky model, with 
eddy-coefficient Cs = 0.17 prevents the flow from developing beyond 
the transitional stage due to excessive dissipation in the early stages of 
the evolution. Finally, the dynamic mixed models were all shown to 
perform about equally well and provide accurate predictions. 

4. LES-ck of a mixing layer 

In this section we will consider LES using the LES-a model. Above, 
in section 2.2, we introduced this model and identified three distinct 
contributions; in fact, the LES-a model contains the explicitly filtered 
nonlinear gradient model (mfj ), the Leray model (mfj) and the com- 
plete LES-a model (mf ). These are defined as 



First we will consider reference LES using these models and compare 
predictions with those obtained with dynamic subgrid models (subsec- 
tion 4.1). Then we focus our attention on the resolved kinetic energy dy- 
namics in subsection 4.2. Finally, in subsection 4.3 we consider (nearly) 
grid-independent LES-a predictions which arise when refining the grid 
while keeping A constant. 

4.1. Reference LES of the mixing layer 

In order to create a point of reference, we consider LES defined on a 
resolution of 32 3 grid-points. This choice represents a significant saving 
compared to the full DNS and places a considerable importance on the 
subgrid fluxes. This resolution was used previously in a comparative 
study of subgrid models in [7] . 

The simulations will be illustrated by considering the evolution of the 
resolved kinetic energy E(t), defined in (67). In addition, we consider 
the momentum thickness S(t), based on filtered variables which quanti- 
fies the spreading of the mean velocity profile. We also investigate the 
Reynolds-stress profiles {w±W2) defined with respect to the fluctuation 
Wi = Vi — (vi). Finally, we incorporate the streamwise kinetic energy 
spectrum in the turbulent regime at t = 80. In this way a number of es- 
sentially different quantities (mean, local, plane averaged) are included 




(68) 
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in the comparisons in order to assess various aspects of the quality of 
the models. 

For all simulations we will use a LES-filter-width A = L/16. On the 
32 3 grid this implies that A/h = 2, i.e., two grid-intervals cover the 
filter-width. Moreover, unless explicitly stated otherwise, the explicit 
filter used in the definition of the LES-a subgrid models will have the 
same width as the LES-filter, i.e., k = 1. The filtering is done using the 
top-hat filter and we adopt the trapezoidal rule to perform the numerical 
integrations. The simulations that will be presented in the following 
subsections correspond to a slightly different initial condition than used 
in section 3.2. The differences are fairly small, but still prevent a direct 
comparison with the filtered DNS results presented in section 3.2. 

Explicit filtering is essential. The proposed subgrid models in 
the a framework each contain the nonlinear gradient model and also 
involve an explicit filtering. As analyzed in section 2.4, the nonlinear 
gradient model, without explicit filtering gives rise to instabilities. These 
instabilities manifest themselves, e.g., by an increase in the resolved 
kinetic energy, instead of the monotonous decrease that is characteristic 
of this relaxing shear layer, cf. figure 3. 




Figure 4- Evolution of resolved kinetic energy for the nonlinear gradient model 
(solid) and the filtered nonlinear gradient model, using n = 1 (dashed) and k = 
2 (dash-dotted) (a). In (b) we show the corresponding results obtained with the 
unfiltered (solid) and filtered full LES-a model. These instabilities are expected on 
grounds discussed in subsection 2.4. 

The question arises whether the explicit filtering can stabilize the sim- 
ulations on this reference grid. In figure 4 we compiled predictions for 
the kinetic energy, obtained with the nonlinear gradient and the full 
LES-a model, both without and with explicit filtering, at different val- 
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ues of the ratio k. We notice that the explicit filtering is essential in 
order to maintain stability of the simulation. It appears that the un- 
filtered LES-a model is even slightly more unstable than the unfiltered 
nonlinear gradient model. We also considered these models at a higher 
resolution of 64 3 grid-points. Consistent with the analysis in section 2.4 
the instability becomes stronger if the grid is refined while keeping the 
LES filter-width A constant. It is seen that the value of k, which defines 
the width of the explicit filter relative to the width of the LES filter, has 
a comparably small effect on the predictions of the nonlinear gradient 
model. The instabilities which arise when using the full LES-a model, 
appear somewhat stronger and, e.g., E even increases in the turbulent 
regime, despite the explicit filtering. This indicates a marginally unsta- 
ble simulation, and the situation improves when k is increased. 

Reference LES-ct predictions. Some basic predictions obtained 
using the three LES-a models will be presented next. These predictions 
will contain errors because of shortcomings in the subgrid parameteriza- 
tions and the numerical treatment. These aspects will be focused upon 
in the next two subsections respectively; here it is our aim to provide an 
impression of the predictions under numerical conditions that are fairly 
common in present-day LES, e.g. A/h = 2. 
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Figure 5. Evolution of the resolved kinetic energy E comparing the following mod- 
els: LES-a (dashed), Leray (solid), filtered nonlinear gradient model (dash-dotted), 
dynamic mixed (+), dynamic eddy- viscosity (o) and no-model (o). We used k = 1. 

In figure 5 we compare the evolution of the resolved turbulent ki- 
netic energy E for a number of subgrid models. We included not only 
predictions corresponding to the three LES-a models, but also the dy- 
namic mixed model, the dynamic eddy-viscosity model and the simula- 
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tion without any subgrid model at all. The subgrid models provide a 
significant improvement compared to the case without a model. From 
previous simulations we know that a fairly close agreement exists be- 
tween filtered DNS data and the dynamic models, as shown in figure 3 
(see [7] for more details). Using the dynamic predictions as point of 
reference here as well, we notice that the Leray and the filtered non- 
linear gradient model provide more accurate predictions than the full 
LES-a model. We also considered the Bardina model and observed that 
the predictions are virtually identical to those obtained with the filtered 
nonlinear gradient model. The Smagorinsky model at Cs = 0.17 was 
used as well and showed too strong dissipation. 
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Figure 6. Evolution of the resolved momentum thickness 8 comparing the follow- 
ing models: LES-a (dashed), Leray (solid), filtered nonlinear gradient model (dash- 
dotted), dynamic mixed (+), dynamic eddy- viscosity (o). We used n = 1. 

The momentum thickness 5 is shown in figure 6. The prediction of 
5 from the full LES-a model is much higher than those obtained with 
the other subgrid models and compared to the dynamic model predic- 
tions as point of reference, it appears too high. The predictions of the 
Bardina similarity model again coincide with the filtered nonlinear gra- 
dient model, and these predictions are somewhat larger than arise from 
the Leray model. All the LES-a models predict S larger than the dy- 
namic models. Since the dynamic predictions slightly underestimate 5 
according to [7] , it appears that the Leray model and the filtered nonlin- 
ear gradient model predict 5 more accurately, compared to filtered DNS 
results, than the other models. 

In figure 7 we collected the Reynolds stress — (ti^ii^)- We observe that 
all three LES-a models predict a considerably higher level of fluctuations 
compared to the dynamic models. The full LES-a model predicts levels 
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Figure 7. Comparison of the Reynolds stress —{1V1W2) at t = 70: LES-a (dashed), 
Leray (solid), filtered nonlinear gradient model (dash-dotted), dynamic mixed (+), 
dynamic eddy- viscosity (o) and no- model (o) . We used k = 1 . 

of fluctuation close to those obtained from the simulation without any 
subgrid model, suggesting that this model introduces too many small 
scales into the solution. Likewise, the filtered gradient model generates 
high levels of fluctuations, while the Leray model is much closer to the 
levels of fluctuation that are found using the dynamic models. 
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Figure 8. Comparison of the streamwise energy spectrum A(k) at t = 80: 
LES-a (dashed), Leray (solid), filtered nonlinear gradient model (dash-dotted), dy- 
namic mixed (+), dynamic eddy- viscosity (o) and no-model (o). We used k = 1. 

We consider the streamwise kinetic energy spectrum in the turbulent 
regime at t = 80, in figure 8. We observe a clear separation of the predic- 
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tions in two groups. The two dynamic models show a strong reduction 
of the smaller scales. In contrast, the full LES-a model displays a spec- 
trum that is quite close to the spectrum of the simulation without any 
subgrid model. This situation improves significantly for the filtered non- 
linear gradient model and finally, the Leray model provides the largest 
attenuation of the small scales among the three LES-a models. 

In summary, the simulations suggest that the full LES-a model does 
not sufficiently reduce the resolved kinetic energy, leads to too large 
momentum-thickness and too high levels of fluctuation, which is ap- 
parent in the spectrum at small scales and snapshots of the solution. 
The filtered nonlinear gradient model performs better than the full 
LES-a model but also over-predicts the smaller scales. In contrast to 
these two models, the Leray model, appears to predict the energy de- 
cay properly, shows accurate momentum-thicknesses and apparently re- 
liable levels of turbulence intensities, as shown also in the spectrum and 
in snapshots of the solution. In order to better understand these pre- 
dictions we turn to the resolved kinetic energy dynamics in the next 
subsection and consider the contribution of the individual terms in the 
models. 

4.2. Resolved kinetic energy dynamics 

In this section we consider the evolution of the resolved kinetic energy 
and determine the type and magnitude of the various subgrid contribu- 
tions. The evolution of E is governed by 




— UidjTij} dx 



^~ 2Re aij a ' ij + Tij9jUi ^ dx ( 71 ) 

where use was made of the identity WijdjUi = ^aJJ ajj. The predicted 
kinetic energy evolution, corresponding to a given LES model, emerges 
by replacing the turbulent stress tensor by its subgrid scale model. We 
notice that the dynamics of E is governed by a purely dissipative term 
arising from the molecular dissipation and a term that is associated 
with the subgrid model. We will consider the resolved energy dynamics 
both for the coarse reference grid of 32 3 grid points and a much finer 
simulation in which we use 96 3 . The latter simulations use the same 
filter-width A = L/16 but correspond to a much higher subgrid reso- 
lution A = GJiles- I n this wa Y we can clarify some of the dynamics 
observed on the coarse grid as well as obtain an impression of the actual 
dynamical consequences associated with the subgrid model. 
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Figure 9. Resolved kinetic energy rate contributions as a function of time: 
LES-a (dashed line dtE m , +: dtE v ), Leray (solid line: dtE m , o: dtE v ), filtered 
nonlinear gradient model (dash-dotted line dtE m , o: dtE v ). We used k = 1 and show 
the results for the 32 3 grid in (a) and for the 96 3 grid in (b) . 



In figure 9 we show the total viscous and subgrid contributions to 
dtE, denoted dfE v and dtE m , respectively. We notice that on the 32 3 
grid the viscous contribution corresponding to the Leray model is quite 
constant in the turbulent regime and the subgrid contribution gradu- 
ally becomes of the same order of magnitude. For the filtered nonlinear 
gradient model we observe a proper dissipation of energy, but slightly 
less than the Leray model. The corresponding viscous flux contribu- 
tion increases considerably in the turbulent regime. Finally, for the full 
LES-a model we observe that the subgrid contribution not only becomes 
less important in the turbulent regime but even changes sign. This can 
readily be associated with the overestimated small scale contributions in 
the solution, as shown in the previous subsection. For the better resolved 
LES the results of the Leray model and the filtered nonlinear gradient 
model are quite comparable and appear more predictable. Moreover, 
all subgrid fluxes are seen to settle and oscillate around some nonzero 
values, indicating perhaps a more regular self-similar development of the 
mixing layer in the turbulent regime. The full LES-a model was found 
to become unstable around t = 70, at this high subgrid resolution. Ap- 
parently, the explicit filtering, which was found to be essential in the 
previous section, in order to stabilize the simulation on the coarse grid, 
is not damping sufficiently well to maintain stability of the LES-a model 
at increased subgrid resolution. 

To further analyse the dynamical behavior, we can look at splitting 
the subgrid contribution into a positive, i.e., forward scatter or dissi- 
pative, contribution and a negative, i.e., backward scatter or reactive 
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Figure 10. Resolved kinetic energy rate contributions: LES-a (dashed line Pj, +: 
Pb), Leray (solid line: Pf,o: Pt ), filtered nonlinear gradient model (dash-dotted line 
Pf,o: Pt)- We used k = 1 and show the results for the 32 3 grid (a) and the 96 3 grid 
(b). 



contribution. To formalize this splitting, we introduce 

W) = I ' \(f + \f\)dx, P f (f) = f hf-\f\) dx (72) 
Jn 2 Jq 2 

to measure the amount of back-scatter and forward scatter of energy 
(Pf) associated with a term represented by /. In figure 10 we collected 
the forward and backward scatter contributions for the LES-a models. 
We observe that all these models predict both forward and backward 
scatter of energy, which sets them apart from simple eddy- viscosity mod- 
els that only provide forward scatter. On the coarse grid (32 3 ) the Leray 
model and the filtered nonlinear gradient model compare fairly well. The 
full LES-a model, however, shows a large amount of back-scatter in the 
turbulent regime and a likewise increased importance of forward scatter. 
On the finer grid the Leray and filtered nonlinear gradient model show a 
balance between forward and backward scatter in the turbulent regime. 

A third decomposition of the total contribution arises in terms of the 
individual subgrid-terms. If we consider, e.g., the full LES-a model, 
written as mfj = Aij + Bij — Cij we may write 

d t E = 3 t E v + 8 t E A + d t E B - d t E c (73) 

with individual contributions due to the viscous fluxes, and the A—B — C 
terms respectively; 

d t E v = ( -TT^-oTJoTj dx , d t E A = f A^djlli dx (74) 
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Figure 11. Resolved kinetic energy rate contributions: full LES-a shown in (a) 
(solid A-term, dashed B-term, dash-dotted C-term), Leray model shown in (b): (solid 
A-term, dashed B-term) We used n = 1 and show the results for the 32 3 grid. 

and similarly for the other terms. In figure II we collected the detailed 
energy-dynamics decomposition corresponding to the two terms which 
make up the Leray model and the three terms that constitute the full 
LES-a model. Notice that figure 9 already contains the single contri- 
bution of the filtered nonlinear gradient model. The Leray model is 
seen to be composed of two terms that both dissipate energy. The full 
LES-a model behaves less regular and we observe that the dissipative B 
contribution is nearly canceled by the reactive C contribution. 

From this analysis of the resolved energy dynamics it seems that the 
Leray model and the filtered nonlinear gradient model provide more 
accurate results and the internal functioning of these models increases 
the robustness of the model. We also applied the Leray model to a flow 
at a ten times higher Reynolds number. Although these latter results are 
still preliminary, it seems that the Leray model provides reliable results, 
even in such very turbulent flows. Further analysis of this regime is 
needed though and this will be published elsewhere [37]. 

In the next section we use the well resolved LES predictions in com- 
bination with the coarse grid simulation results to assess the influence 
of the spatial discretization scheme on the evolution of the flow. 

4.3. Toward grid-independent LES-a 

The reference simulations considered above, are executed on a fairly 
coarse grid which corresponds to a ratio between filter-width A and 
grid-spacing h of A/h = 2. In many present-day LES, even a ratio of 
A/h = 1 is frequently used. These choices usually arise from considera- 
tions of available computational resources, but at the same time imply 
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that the smallest resolved scales of size on the order of A are not accu- 
rately represented in the numerical treatment. Hence, there is a strong 
possibility that the marginal subgrid resolution influences the dynamical 
properties of the simulated flow. 

In order to assess this discretization effect, we will compare the refer- 
ence LES with simulations at a much higher subgrid resolution. In this 
way, the effect of subgrid modeling is better represented numerically, 
while it remains of the same magnitude as in the coarse grid simulation. 
This allows to isolate the dynamic effects of the spatial discretization in 
the modeled equations. 




Figure 12. Resolved kinetic energy using the Leray model (a) and the filtered 
nonlinear gradient model (b): solid line (4th order method, resolution 96 3 ), dashed 
line (4th order method, resolution 64 ), dash-dotted line (4th order method, resolution 
32 3 ) and dotted-line (2nd order method, resolution 32 3 ). 

We compare simulations on 32 3 , 64 3 and 96 3 grid-points and focus 
our attention on the Leray model and the filtered nonlinear gradient 
model. In figure 12 we compare the predicted resolved kinetic energy 
obtained with the second and the fourth order accurate spatial discretiza- 
tion method. The subgrid resolution corresponding to these three grids 
is A/h = 2, 4 and 6 respectively. We observe a very close agreement 
between the predictions using the fourth order accurate method and 
A/h = 4 and 6. This suggests that a mean flow quantity such as the 
resolved kinetic energy is well represented using A/h = 4. Moreover, we 
notice that on the coarsest grid, the accuracy of the prediction based 
on the second order method compares closely to that obtained with the 
fourth order method. Apparently, if the dynamic effects of the spatial 
discretization errors are quite large, a lower order method can be com- 
petitive with a higher order method. For both models the reliability of 
the predictions on the coarse grid are affected considerably by the coarse- 
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ness of the subgrid resolution. In both situations, and for both spatial 
discretizations the effect of the discretization error is seen to enhance 
the reduction of the resolved kinetic energy. 




Figure 13. Resolved momentum-thickness using the Leray model (a) and the filtered 
nonlinear gradient model (b): solid line (4th order method, resolution 96 3 ), dashed 
line (4th order method, resolution 64 3 ), dash-dotted line (4th order method, resolution 
32 3 ) and dotted-line (2nd order method, resolution 32 3 ). 



To further establish the convergence, we show the momentum-thickness 
in figure 13. We observe that the convergence is clear for the Leray 
model and that a value of A/h = 4 corresponds to reliable predictions 
for both subgrid models considered, although the sensitivity of the mo- 
mentum thickness is larger than that of the resolved kinetic energy. Re- 
garding the results for the best resolved simulations, we observe that 
the momentum-thickness develops very nearly linearly with time in the 
Leray model, while a slight reduction of the growth-rate predicted by 
the nonlinear gradient model is seen in the turbulent regime. 

Finally, we show the spectra obtained with the Leray and the fil- 
tered nonlinear gradient model on the selected grids in figure 14. We 
notice a general resemblance between the results obtained with both 
models. As the subgrid resolution is increased, a larger portion of the 
spectrum is better resolved, cf. the spectra obtained on 64 3 and 96 3 
grid-points. Moreover, the differences due to the use of the second order 
or the fourth order accurate methods are expressed very clearly in the 
spectra; a strong reduction in energy in the higher wavenumbers on the 
32 3 grid results when using the second order method. This is consistent 
with the stronger attenuation of the high wavenumbers in the second 
order method, compared to the fourth order accurate method. 
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Figure 14- Streamwise kinetic energy spectrum using the Leray model (a) and the 
filtered nonlinear gradient model (b): solid line (4th order method, resolution 96 3 ), 
dashed line (4th order method, resolution 64 3 ), dash-dotted line (4th order method, 
resolution 32 3 ) and marker 'o' (2nd order method, resolution 32 3 ). 

5. Concluding remarks 

In the a-framework we derived a new subgrid closure for the turbulent 
stress tensor in LES, by using the Kelvin theorem applied to a filtered 
transport velocity. The proposed LES-a subgrid model was shown to 
contain two other subgrid parameterizations, i.e., the nonlinear gradient 
model and a model corresponding to the Leray regularization of Navier- 
Stokes dynamics. Moreover, the LES-a model stress tensor contains an 
explicit filtering in its definition, which sets it apart from other subgrid 
models in literature. It was shown that this explicit filtering is essential 
for the LES-a models; without it, the simulations develop a finite time 
instability. This instability was also observed in an analysis of the viscous 
one-dimensional Burgers equation and appears to be associated with the 
nonlinear gradient term. 

The flow in a turbulent mixing layer was considered, in order to test 
the capabilities of the three 'nested' LES-a models. Through a compar- 
ison with dynamic (mixed) models, we inferred that the Leray model 
provides particularly accurate predictions. The filtered nonlinear gradi- 
ent model in turn, compares well with the Leray model and corresponds 
closely to the Bardina similarity model for the flows considered. The full 
LES-a model was seen to generate too many small scales in the solution 
and correspondingly poorer predictions, e.g., too large growth-rate, too 
high levels of turbulence intensities, etc. An analysis of the resolved 
kinetic energy dynamics showed that the full LES-a model contains two 
competing contributions which may tend to cancel and, thus, destabilize 
some simulations involving this model. In particular, this tendency im- 
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plied that simulations using the full LES-a model were unstable in our 
simulations of turbulent mixing at high subgrid resolution. 

Apart from accuracy and a certain degree of numerical robustness, the 
computational overhead associated with a subgrid model is an element of 
importance in evaluating simulation methods. The computational effort 
associated with the three LES-a models is considerably lower than that 
of dynamic (mixed) models. This is primarily a result of the reduction 
in the number of explicit filtering operations required to evaluate the 
LES-a model. Moreover, the accuracy of the predictions is higher for 
the Leray model than for any of the other subgrid models considered. 
Consequently, the Leray model is favored in this study and holds promise 
for applications to even more complex and demanding flow problems. 
Preliminary results at significantly higher Reynolds number suggest that 
the Leray model performs well also in this case. 

Finally, we also considered the contribution to the dynamics arising 
from the spatial discretization method at coarse subgrid resolution. In 
general, the role of numerical methods in relation to LES has not yet 
been sufficiently clarified to determine unambiguously whether the accu- 
racy of predictions is restricted because of shortcomings in the subgrid 
model, or whether this inaccuracy is due to spatial discretization ef- 
fects. Resolving this ambiguity and determining the main sources of 
error would help in finding the best strategy for employing computa- 
tional resources in an LES. In one strategy, a grid- independent solution 
of the modeled LES equations at fixed filter-width is sought and one 
only assigns computational resources for reducing numerical errors by 
increasing the subgrid resolution ratio A/h. This approach was applied 
here and used to evaluate the additional dissipation that arises from the 
spatial finite volume discretizations at either second order, or fourth or- 
der accuracy. This additional dissipation is associated with the implicit 
filtering effect of the small flow features represented on the grid. 
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